LP és MIP feladatok megoldása Pythonnal

Lineáris és egészértékű programozás

Lineáris programozás: Lineáris függvény optimalizása lineáris feltételek mellett. Legyen $A\in\mathbb{R}^{n\times n}$ feltételi mátrix, $b\in\mathbb{R}^m$ korlátozóvektor és $c\in\mathbb{R}^n$ célfüggvény. Olyan $x\in\mathbb{R}^n$ vektort keresünk, ami a feltételrendszert kielégítő vektorok közül a célfüggvény szerint maximális. Az $x_i$ változókat döntési változónak hívjuk. A primál lineáris program (LP) az alábbi:

\begin{equation} \begin{array}{rrcl} \max & cx\\ \mathrm{s.t.} &Ax &\leq &b\\ &x&\geq&0. \end{array} \end{equation}

Erre gondolhatunk a következő formában is:

\begin{equation} \begin{array}{rrcll} \max & \sum_{j=1}^nc_jx_j\\ \mathrm{s.t.} &\sum_{j=1}^na_{ij}x_j &\leq &b_i&\forall i=1,\dots,m\\ &x_j&\geq&0&\forall j=1,\dots,n. \end{array} \end{equation}

Ekkor az egyenlőtlenségek (hiper-)féltereket jelentenek, ezek metszete az a halmaz, amely felett keressük a célfüggvény optimumát. Ezt primál poliédernek nevezzük.A fenti program duálisa:

\begin{equation} \begin{array}{rrcl} \min & yb\\ \mathrm{s.t.} &yA &\geq &c\\ &y&\geq&0, \end{array} \end{equation}

ahol $y\in\mathbb{R}^m$. Ez is átírható szummás alakba: \begin{equation} \begin{array}{rrcll} \min & \sum_{i=1}^m y_ib_i\\ \mathrm{s.t.} &\sum_{i=1}^m y_ia_{ij} &\geq &c_j& \forall j=1,\dots,n\\ &y_i&\geq&0&\forall i=1,\dots,m. \end{array} \end{equation}

Az lineáris programozásra vonatkozó dualitás tétel alapján ha a primál és duál közül az egyiknek van véges optimuma, akkor a másiknak is van ($x^*,y^*$), és ekkor

$$ y^*b=cx^*. $$

Az lineáris programozás geometriai értelmezésében tetszőleges lineáris célfüggvény esetén létezik olyan optimális megoldás, amely a poliéder csúcsa. Ez az alapja a Simplex módszernek.

Egészértékű programozás: Az egészértékű programozás alatt a lineáris programozásnak azt a változatát értjük, amelyben néhány döntési változóra egészértékűséget követelünk meg. Például legyen $I\subseteq 1,\dots,n$, ekkor a következő optimalizálási feladat egy (vegyes) egészértékű program (IP vagy MIP):

\begin{equation} \begin{array}{rrcll} \max & \sum_{j=1}^nc_jx_j\\ \mathrm{s.t.} &\sum_{j=1}^na_{ij}x_j &\leq &b_i&\forall i=1,\dots,m\\ &x_j&\geq&0&\forall j=1,\dots,n\\ &x_j&\in&\mathbb{Z}&\forall j\in I. \end{array} \end{equation}

Algoritmusok:

Lineáris programozás:

  • Simplex módszer (primál, duál, különböző pivotálási szabályok)
  • Belsőpontos módszerek:
    • Ellipsoid módszer
    • Affin skálázás
    • Barrier módszerek

Egészértékű programozás:

  • Branch-and-Bound
  • Vágósíkos eljárások
  • Branch-and-Cut
  • Branch-and-Price

Solverek

Open source Kereskedelmi
GLPK CPLEX
CLP Gurobi
CBC FICO Xpress
stb. stb.

A legtöbb kereskedelmi solverhez elérhető ingyenes academic license.

Solverek használata Pythonban

Pythonban solverek használatára négy példát fogunk nézni:

  • SciPy LP solverét
  • A solver API-ján keresztül
  • PyOmo-val
  • PuLP segítségével

SciPy LP solver

A scipy package optimize moduljából a linprog függvényt használjuk. Mikor jó ez? Ha valaki olyan alakú feladatot akar megoldani, hogy

\begin{equation} \begin{array}{rrcl} \max & cx\\ \mathrm{s.t.} &A_{\leq}x &\leq &b_{\leq}\\ &A_=x&=&b_=\\ &x&\geq&l \\ &x&\leq&u, \end{array} \end{equation}

ahol az $A_\leq,A_=$ feltételi mátrixok, $b_\leq,b_=$ korlátozó vektorok, $l,u$ alsó és felső korlátok és $c$ költségfüggvény explicit adottak.

A függvény bemenetként numpy arrayeket vár, és kimenetként egy dictionaryt ad vissza, amely a megoldásról tartalmaz adatokat ($x^*$ optimális megoldást, optimum értéket, iterációk számát stb.). A függvény hívásakor lehet a megoldási módszert is kiválasztani, itt lehet válogatni különböző simplex variácók közül illetve belső pontos módszer használata a default megoldási eljárás. (Az implementáció a HiGHS solveré.)

Ez sajnos nem tud IP-t megoldani, azonban ennek a segítségével lehetséges például egy LP alapú Branch-and-Bound módszert készíteni.

A függvény szignatúrája:

scipy.optimize.linprog(c,
                       A_ub=None,
                       b_ub=None,
                       A_eq=None,
                       b_eq=None,
                       bounds=None,
                       method='interior-point',
                       callback=None,
                       options=None,
                       x0=None)

Ehhez útmutató: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html

Keverési feladat

A feladat: Példaként egy keverési feladaton fogjuk megnézni, hogy hogyan kellene a SciPy LP solverével modellezni és megoldani. Tegyük fel, hogy egy macskatápot szeretnénk kikeverni.

Alapanyag Fehérje (g) Zsír (g) Rost (g) Só (g)
Csirke (g) 0.1 0.08 0.001 0.002
Marha (g) 0.2 0.1 0.005 0.005
Birka (g) 0.15 0.11 0.003 0.007
Rizs (g) 0.000 0.01 0.1 0.002
Korpa (g) 0.04 0.01 0.15 0.008
Zselatin (g) 0 0 0 0

Az alapanyagok költsége:

  Csirke (g) Marha (g) Birka (g) Rizs (g) Korpa (g) Zselatin (g)
Költség ($) 0.013 0.008 0.010 0.002 0.005 0.001

Tudjuk továbbá, hogy egy konzerv macskatáp pontosan 100 gram, és a tápanyagokból megfelelő mennyiséget tartalmaz, amit az alábbi táblázatban láthatunk:

Tápanyag Mennyiség (g)
Fehérje min. 8
Zsír min. 6
Rost max. 2
max. 0.4

A cél elkészíteni egy keveréket az alapanyagokból, amely megfelel a fenti tápanyagkövetelményeknek, és a lehető legolcsóbb.

A modell: Legyenek a döntési változóink $x_1,x_2,x_3,x_4,x_5,x_6$, ahol a változók jelentése az alábbi:

Változó Alapanyag (g)
$x_1$ Csirke
$x_2$ Marha
$x_3$ Birka
$x_4$ Rizs
$x_5$ Korpa
$x_6$ Zselatin

Célunk a lehető legolcsóbban előállítani a macskatápot, ezért a célfüggvény

$$ \min 0.013 x_1+0.008 x_2 + 0.01 x_3+0.002 x_4+0.005 x_5+0.001 x_6. $$

Tudjuk továbbá, hogy pontosan 100 g a nettó tömege egy konzerv tápnak, ezért azt az egyenletet felvehetjük a feltételeink listájába, hogy

$$ x_1+x_2+x_3+x_4+x_5+x_6 = 100. $$

Az alapanyagok tápanyagtartalma és a keverékre vonatkozó tápanyagkövetelmények alapján azt kapjuk, hogy a fehérjetartalmora vonatkozó korlát

$$ 0.1x_1+0.2x_2+0.15x_3+0.04x_5 \geq 8, $$

ahol a rizst és a zselatint kihagytam a korlátból, mivel a fehérjetartalmuk 0. A zsírtartalomra vonatkozó korlát

$$ 0.08x_1+0.1x_2+0.11x_3+0.01x_4+0.01x_5\geq 6, $$

a rosttartalomra vonatkozó korlát

$$ 0.001x_1+0.005x_2+0.003x_3+0.1x_4+0.15x_5\leq 2, $$

a sótartalomra vonatkozó korlát pedig

$$ 0.002x_1+0.005x_2+0.007x_3+0.002x_4+0.008x_5 \leq 0.4. $$

Átírás mátrix alakra: A scipy.optimize.linprog által támogatott mátrixos alakra átírjuk a fenti korlátokat. A célfüggvény:

In [41]:
import numpy as np
from scipy.optimize import linprog
In [42]:
c = np.array([0.013, 0.008, 0.010, 0.002, 0.005, 0.001])

Az $A_=$ feltételi mátrix és a $b_=$ korlátozóvektor:

In [43]:
Aeq = np.array([[1,1,1,1,1,1]])
beq = np.array([100])

Az $A_\leq$ feltételi mátrix létrehozásánál figyelnünk kell arra, hogy a függvény $A_\leq x\leq b_\leq$ alakban várja a feltételeket, ezért azokat a sorokat, amelyekben $\geq$ szerepel, meg kell szoroznunk $-1$-gyel:

In [69]:
Aieq = np.array([[-0.100, -0.200, -0.150,  0.000, -0.040,  0.000],
                 [-0.080, -0.100, -0.110, -0.010, -0.010,  0.000],
                 [ 0.001,  0.005,  0.003,  0.100,  0.150,  0.000],
                 [ 0.002,  0.005,  0.007,  0.002,  0.008,  0.000]])

bieq = np.array([-8.0,
                 -6.0,
                  2.0,
                  0.4])

Ezután nincs más dolgunk, mint odaadni a függvénynek mint input:

In [70]:
blend = linprog(c, Aieq, bieq, Aeq, beq)
blend
Out[70]:
     con: array([6.29560759e-09])
     fun: 0.5200000026527398
 message: 'Optimization terminated successfully.'
     nit: 8
   slack: array([3.99999998e+00, 8.83075391e-09, 1.69999997e+00, 9.99999972e-02])
  status: 0
 success: True
       x: array([8.63787793e-08, 5.99999994e+01, 5.64719718e-07, 2.01189816e-08,
       2.26632479e-07, 3.99999997e+01])

Ez azt jelenti, hogy nagyjából 60 gram marhahús és 40 gram zselatin felhasználásával lehet a legolcsóbban, \$0.52 dollárért előállítani egy konzerv macskatápot.

Solver API (FICO Xpress)

A legtöbb kereskedelmi és open source solvernek a funkcionalitását el tudjuk érni a népszerűbb programozási nyelveken (mint például C, C++, Java, Python stb.) keresztül is, úgynevezett API (Application Programming Interface). Ehhez szükséges általában, hogy telepítve legyen a solver és a solver és python kommunikációjához szükséges python csomag az adott számítógépre.

Az előnye ennek, hogy maximálisan ki lehet használni a solver nyújtotta funkcionalitást, a hátránya pedig hogy ismerni kell hozzá a solver API-ját (ami solverenként eltérő lehet, kisebb-nagyobb hasonlóságokkal).

Példaként a FICO Xpress solvert fogjuk megnézni (de nézhetnénk akár Gurobit, Cplexet vagy GLPK-t is). Ehhez útmutató: https://www.msi-jp.com/xpress/learning/square/01-python-interface.pdf

Xpressben a modellt egy problem objektumban tároljuk, ami döntési változók és korlátok sokasága. Itt már nem kell törődnünk azzal, hogy megfelelő alakú mátrixként adjuk a solvernek a korlátjainkat, a korlátok mátrix-szá alakítását már a program végzi el helyettünk, ezzel sokkal kényelmesebb és rugalmasabb modellezést téve lehetővé. (Természetesen elfogad mátrix formában is feltételeket.)

Lássuk az előző feladatot Xpress-szel megoldva!

In [46]:
# Adatok
c = np.array([0.013, 0.008, 0.010, 0.002, 0.005, 0.001])

Aeq = np.array([[1,1,1,1,1,1]])

beq = np.array([100])

Aieq = np.array([[-0.1, -0.2, -0.15, 0, -0.04, 0],
                 [-0.08, -0.1 , -0.11, -0.01, -0.01, 0],
                 [0.001, 0.005, 0.003, 0.1, 0.15, 0],
                 [0.002, 0.005, 0.007, 0.002, 0.008, 0]])

bieq = np.array([-8,
                 -6,
                  2,
                  0.4])
In [71]:
import xpress as xp

# Létrehozom a probléma példányt
prob = xp.problem()

# Létrehozok változókat
x = np.array([xp.var() for _ in range(6)])

# Hozzáadom a problémához a változókat 
prob.addVariable(x)

# Beállítom a célfüggvényt a költség minimalizálására
prob.setObjective(xp.Dot(c, x), sense=xp.minimize)

# Hozzáadom a nettó tömegre vonatkozó egyenletet mint korlát
prob.addConstraint(xp.Dot(Aeq, x) == beq)

# Hozzáadom a tápanyagtartalomra vonatkozó korlátokat
prob.addConstraint(xp.Dot(Aieq, x) <= bieq)

# Megoldom
prob.solve()

# Visszatérek a megoldással
prob.getSolution(x)
FICO Xpress v8.11.3, Hyper, solve started 12:40:58, Oct 7, 2021
Heap usage: 339KB (peak 339KB, 604KB system)
Minimizing LP noname with these control settings:
OUTPUTLOG = 1
Original problem has:
         5 rows            6 cols           25 elements
Presolved problem has:
         5 rows            6 cols           25 elements
Presolve finished in 0 seconds
Heap usage: 340KB (peak 352KB, 606KB system)

Coefficient range                    original                 solved        
  Coefficients   [min,max] : [ 1.00e-03,  1.00e+00] / [ 4.00e-03,  1.00e+00]
  RHS and bounds [min,max] : [ 4.00e-01,  1.00e+02] / [ 8.00e+00,  1.00e+02]
  Objective      [min,max] : [ 1.00e-03,  1.30e-02] / [ 1.00e-03,  1.30e-02]
Autoscaling applied standard scaling

 
   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time
     0           .000000      D      3     0        .000000     0
     2           .520000      D      0     0        .000000     0
Uncrunching matrix
Optimal solution found
Dual solved problem
  2 simplex iterations in 0s

Final objective                       : 5.200000000000000e-01
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max dual violation        (abs/rel) :       0.0 /       0.0
  Max complementarity viol. (abs/rel) :       0.0 /       0.0
Out[71]:
array([ 0., 60.,  0.,  0.,  0., 40.])

Mi történne, ha hozzá akarnék adni még egy olyan korlátot, hogy a marhahús és zselatin össztömege 100 gramban nem haladhatja meg a 90 gramot? Módosítanom kellene a kezdeti feltételi mátrixot, amihez már nem akarok nyúlni. Szerencsére hozzáadhatok még egy korlátot könnyedén:

In [72]:
# Hozzáadok egy új korlátot
prob.addConstraint(x[1]+x[5] <= 90)

# És egoldom újra
prob.solve()
prob.getSolution(x)
FICO Xpress v8.11.3, Hyper, solve started 12:43:40, Oct 7, 2021
Heap usage: 2610KB (peak 2610KB, 606KB system)
Minimizing LP noname with these control settings:
OUTPUTLOG = 1
Original problem has:
         6 rows            6 cols           27 elements
Presolved problem has:
         6 rows            6 cols           27 elements
Presolve finished in 0 seconds
Heap usage: 2610KB (peak 2610KB, 606KB system)

Coefficient range                    original                 solved        
  Coefficients   [min,max] : [ 1.00e-03,  1.00e+00] / [ 4.00e-03,  1.00e+00]
  RHS and bounds [min,max] : [ 4.00e-01,  1.00e+02] / [ 8.00e+00,  1.00e+02]
  Objective      [min,max] : [ 1.00e-03,  1.30e-02] / [ 1.00e-03,  1.30e-02]
Autoscaling applied standard scaling

 
   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time
     0           .520000      D      1     0        .000000     0
     1           .523000      D      0     0        .000000     0
Uncrunching matrix
Optimal solution found
Dual solved problem
  1 simplex iterations in 0s

Final objective                       : 5.230000000000000e-01
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max dual violation        (abs/rel) :       0.0 /       0.0
  Max complementarity viol. (abs/rel) :       0.0 /       0.0
Out[72]:
array([ 0., 59.,  0., 10.,  0., 31.])

Sőt, a változókat nem kötelező tömbben vagy listában tárolni, bármilyen iterálható objektumban (vagy akár önállóan is) elérhetők. A fenti példát átírom egy sokkal átláthatóbb alakra:

In [73]:
# Létrehozom a változókat és egy dictionaryben eltárolom őket
x = {"CSIRKE"    : xp.var(),
     "MARHA"     : xp.var(),
     "BIRKA"     : xp.var(),
     "RIZS"      : xp.var(),
     "KORPA"     : xp.var(),
     "ZSELATIN"  : xp.var()}

# Létrehozom a problémát
prob = xp.problem(x)

# Beállítom a célfüggvényt
prob.setObjective(0.013*x["CSIRKE"]+0.008*x["MARHA"]+0.1*x["BIRKA"]+0.002*x["RIZS"]+0.002*x["KORPA"]+0.001*x["ZSELATIN"],
                  sense=xp.minimize)

# Hozzáadom a nettó tömegre vonatkozó korlátot
prob.addConstraint(xp.Sum(x) == 100)

# Hozzáadom a tápanyagtartalomra vonatkozó korlátokat
prob.addConstraint(0.1*x["CSIRKE"]+0.2*x["MARHA"]+0.15*x["BIRKA"]+0.04*x["KORPA"] >= 8) # Fehérje

prob.addConstraint(0.08*x["CSIRKE"]+0.1*x["MARHA"]+0.11*x["BIRKA"]+0.01*x["RIZS"]+0.01*x["KORPA"] >= 6) # Zsír

prob.addConstraint(0.001*x["CSIRKE"]+0.005*x["MARHA"]+0.003*x["BIRKA"]+0.1*x["RIZS"]+0.15*x["KORPA"] <= 2) # Rost

prob.addConstraint(0.002*x["CSIRKE"]+0.005*x["MARHA"]+0.007*x["BIRKA"]+0.002*x["RIZS"]+0.008*x["KORPA"] <= 0.4) # Só

# Megoldom
prob.solve()

# Visszaadom a megoldást
prob.getSolution(x)
FICO Xpress v8.11.3, Hyper, solve started 12:49:09, Oct 7, 2021
Heap usage: 340KB (peak 340KB, 606KB system)
Minimizing LP noname with these control settings:
OUTPUTLOG = 1
Original problem has:
         5 rows            6 cols           25 elements
Presolved problem has:
         5 rows            6 cols           25 elements
Presolve finished in 0 seconds
Heap usage: 341KB (peak 353KB, 608KB system)

Coefficient range                    original                 solved        
  Coefficients   [min,max] : [ 1.00e-03,  1.00e+00] / [ 4.00e-03,  1.00e+00]
  RHS and bounds [min,max] : [ 4.00e-01,  1.00e+02] / [ 8.00e+00,  1.00e+02]
  Objective      [min,max] : [ 1.00e-03,  1.00e-01] / [ 1.00e-03,  1.00e-01]
Autoscaling applied standard scaling

 
   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time
     0           .000000      D      3     0        .000000     0
     2           .520000      D      0     0        .000000     0
Uncrunching matrix
Optimal solution found
Dual solved problem
  2 simplex iterations in 0s

Final objective                       : 5.200000000000000e-01
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max dual violation        (abs/rel) :       0.0 /       0.0
  Max complementarity viol. (abs/rel) :       0.0 /       0.0
Out[73]:
{'CSIRKE': 0.0,
 'MARHA': 60.0,
 'BIRKA': 0.0,
 'RIZS': 0.0,
 'KORPA': 0.0,
 'ZSELATIN': 40.0}

Eddig nem sok olyat láttunk, amit ne tudtunk volna kevésbé kényelmesen, de a scipy segítségével hasonlóan megoldani. (Nyilván nagyobb feladatpéldányokban már az LP megoldási ideje is sokkal rövidebb a kereskedelmi solverek esetén.)

A legtöbb solver egészértékű programozási feladatot (IP) is tud kezelni, és különböző branch-elési és vágásgenerálási eljárások alkalmazásával sokkal hatékonyabban, mintha nekünk kellene megírni pl egy LP alapú B&B-t ahol csak az Xpress LP solverét használhatnánk.

Nézzünk egy újabb példát.

Szendvicskészítési feladat

A feladat: Tegyük fel, hogy szendvicseket akarunk készíteni, és tegyük fel továbbá, hogy kétféle szendvicset tudunk csinálni, továbbá a szükséges alapanyagokból véges sok áll csak a rendelkezésünkre (leszámítva mondjuk a kenyeret). Továbbá feltételezzük, hogy csak egész számú szendvics készítésére van módunk.

A szendvicstípusok:

Típus Sajt (szelet) Sonka (szelet) Szalámi (szelet) Uborka (karika) Paradicsok (karika)
1 4 1 5 1 0
2 3 3 2 0 1

Az alapanyagmennyiségek:

Alapanyag Mennyiség Egység
Sajt 141 szelet
Sonka 85 szelet
Szalámi 159 szelet
Uborka 35 karika
Paradicsom 30 karika

Ezen alapanyagokból szeretném a két szendvicstípusból a lehető legtöbbet elkészíteni.

A modell: Legyen $x_1, x_2$ az egyes és kettes típusú szendvicshez tartozó egészértékű döntési változók.

A célfüggvény:

$$ \max x_1+x_2 $$

A korlátok az alapanyagok összmennyiségére vonatkoznak, így minden alapanyaghoz felírható egy egyenlőtlenség. A sajtból összesen 141 szelet van, az egyes típusú szendvics elhasznál szendvicsenként 4-et, a kettes 3-at

$$ 4x_1+3x_2 \leq 141, $$

a sonkából 85 szelet van, az egyes 1-et a kettes 3-at használ

$$ x_1+3x_2 \leq 85, $$

a szalámiból 159 szelet van, az egyes 5-öt használ, a kettes 2-t

$$ 5x_1+2x_2 \leq 159, $$

az uborkából 35 karika van, és csak az egyes használja, szendvicsenként 1-et, míg a paradicsomból 30 karika van, és csak a kettes használja, szendvicsenként 1-et

\begin{equation} \begin{array}{rcl} &x_1 &\leq& 35 \\ &x_2 &\leq& 30. \end{array} \end{equation}

A megoldása:

In [76]:
# Létrehozom a két szendvicshez tartozó változókat
x1 = xp.var(name="1.típus",     # elnevezem a változót
            vartype=xp.integer, # egészértékű legyen a változó
            ub=35,              # a felső korlátja legyen 35
            lb=0)               # az alsó korlátja legyen 0

x2 = xp.var(name="2.típus",
            vartype=xp.integer,
            ub=30,
            lb=0)

# Létrehozok egy problémapéldányt ezzel a két változóval
prob = xp.problem(x1, x2)

# Beállítom a célfüggvényt
prob.setObjective(x1+x2, sense=xp.maximize)

# Hozzáadom a korlátokat
prob.addConstraint(4*x1+3*x2 <= 141) # Sajt
prob.addConstraint(  x1+3*x2 <= 85)  # Sonka
prob.addConstraint(5*x1+2*x2 <= 159) # Szalámi

# Megoldom
prob.solve()

# Visszaadom a megoldást
prob.getSolution(x1, x2)
FICO Xpress v8.11.3, Hyper, solve started 12:56:36, Oct 7, 2021
Heap usage: 338KB (peak 338KB, 630KB system)
Maximizing MILP noname with these control settings:
OUTPUTLOG = 1
Original problem has:
         3 rows            2 cols            6 elements         2 globals
Presolved problem has:
         3 rows            2 cols            6 elements         2 globals
LP relaxation tightened
Presolve finished in 0 seconds
Heap usage: 364KB (peak 367KB, 632KB system)

Coefficient range                    original                 solved        
  Coefficients   [min,max] : [ 1.00e+00,  5.00e+00] / [ 5.00e-01,  1.50e+00]
  RHS and bounds [min,max] : [ 3.00e+01,  1.59e+02] / [ 1.00e+00,  4.25e+01]
  Objective      [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]
Autoscaling applied standard scaling

Will try to keep branch and bound tree memory usage below 8.8GB
 *** Heuristic solution found:     3.000000      Time: 0 ***
 *** Heuristic solution found:    29.000000      Time: 0 ***
Starting concurrent solve with dual

 Concurrent-Solve,   0s
            Dual        
    objective   dual inf
 D  40.777778   .0000000
------- optimal --------
Concurrent statistics:
      Dual: 2 simplex iterations, 0.00s
Optimal solution found
 
   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time
     2         40.777778      D      0     0        .000000     0
Dual solved problem
  2 simplex iterations in 0s

Final objective                       : 4.077777777777777e+01
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max dual violation        (abs/rel) :       0.0 /       0.0
  Max complementarity viol. (abs/rel) :       0.0 /       0.0

Starting root cutting & heuristics
 
 Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time
a           31.000000    40.777778      3                 23.98%       0      0
b           33.000000    40.777778      4                 19.07%       0      0
g           35.000000    40.777778      5                 14.17%       0      0
q           40.000000    40.777778      6                  1.91%       0      0
 *** Search completed ***
Uncrunching matrix
Final MIP objective                   : 4.000000000000000e+01
Final MIP bound                       : 4.000000000000000e+01
  Solution time / primaldual integral :         0s/ 69.532045%
  Number of solutions found / nodes   :         6 /         1
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max integer violation     (abs    ) :       0.0
Out[76]:
(19.0, 21.0)

Pyomo

A Pyomo (Python Optimizing Modeling Objects) egy Python alapú optimalizálási modellező nyelv, lineáris és egészértékű programokon kívül még nagyon sokféle matematikai optimalizálási probléma modellezésére alkalmas.

Honlap: http://www.pyomo.org/

Dokumentáció: https://pyomo.readthedocs.io/en/stable/index.html

Pyomo-ban kétféle modellt lehet készíteni, egyik a konkrét modell (ConcreteModel), másik az absztrakt modell (AbstractModel). A konkrét modell esetében az adatoknak, adatok számának, változók számának már a modell megalkotásakor rendelkezésre kell állnia, míg egy absztrakt modell ezeket paraméterként kezeli.

A Pyomo nagyon hasonlít az algebrai modellező nyelvekre, mint az AMPL, GAMS, Mosel stb. Ez a szintaxison kívül abban nyilvánul meg, hogy a solver nem a Pyomohoz tartozik, hanem a Pyomo egy tőle független, külső (kereskedelmi vagy open source) solvert hív meg a megoldás megtalálásához.

Amik egy modellhez tartozhatnak:

  • Var-ok
  • Set-ek
  • Param-ok
  • Constraint-ek
  • Objective

Nézzük meg az előző keverési feladatot konkrét modellként!

In [77]:
import pyomo.environ as pyo

# Inicializálok egy konkrét modellt
model = pyo.ConcreteModel()

# Hozzáadok egy változót, amit a CSIRKE, MARHA ... stringekkel indexelek
model.x = pyo.Var(["CSIRKE","MARHA","BIRKA","RIZS","KORPA","ZSELATIN"], domain=pyo.NonNegativeReals)

# Hozzáadom a költségfüggvényt
model.cost = pyo.Objective(
    expr = 0.013*model.x["CSIRKE"]+0.008*model.x["MARHA"]+0.1*model.x["BIRKA"]+0.002*model.x["RIZS"]+0.002*model.x["KORPA"]+0.001*model.x["ZSELATIN"])

# Hozzáadom a nettó tömegre vonatkozó korlátot
model.weight = pyo.Constraint(
    expr = model.x["CSIRKE"]+model.x["MARHA"]+model.x["BIRKA"]+model.x["RIZS"]+model.x["KORPA"]+model.x["ZSELATIN"] == 100)

# Hozzáadom a tápanyagtartalomra vonatkozó korlátokat
model.protein = pyo.Constraint(
    expr = 0.1*model.x["CSIRKE"]+0.2*model.x["MARHA"]+0.15*model.x["BIRKA"]+0.04*model.x["KORPA"] >= 8)

model.fat = pyo.Constraint(
    expr = 0.08*model.x["CSIRKE"]+0.1*model.x["MARHA"]+0.11*model.x["BIRKA"]+0.01*model.x["RIZS"]+0.01*model.x["KORPA"] >= 6)

model.fibre = pyo.Constraint(
    expr = 0.001*model.x["CSIRKE"]+0.005*model.x["MARHA"]+0.003*model.x["BIRKA"]+0.1*model.x["RIZS"]+0.15*model.x["KORPA"] <= 2)

model.salt = pyo.Constraint(
    expr = 0.002*model.x["CSIRKE"]+0.005*model.x["MARHA"]+0.007*model.x["BIRKA"]+0.002*model.x["RIZS"]+0.008*model.x["KORPA"] <= 0.4)

# Inicializálok egy solvert, ami az Xpresst fogja használni
solver = pyo.SolverFactory("xpress")

# Megoldom a modellt
solver.solve(model)

# Kiírom az összetételt
print("Költség  :", model.cost(), "$")
print("--------------------")
print("Összetétel")
print("Csirke   :", model.x["CSIRKE"](), "g")
print("Marha    :", model.x["MARHA"](), "g")
print("Birka    :", model.x["BIRKA"](), "g")
print("Rizs     :", model.x["RIZS"](), "g")
print("Korpa    :", model.x["KORPA"](), "g")
print("Zselatin :", model.x["ZSELATIN"](), "g")
print("--------------------")
print("Tápanyagtartalom")
print("Fehérje  :", model.protein(), "g")
print("Zsír     :", model.fat(), "g")
print("Rost     :", model.fibre(), "g")
print("Só       :", model.salt(), "g")
Költség  : 0.52 $
--------------------
Összetétel
Csirke   : 0.0 g
Marha    : 60.0 g
Birka    : 0.0 g
Rizs     : 0.0 g
Korpa    : 0.0 g
Zselatin : 40.0 g
--------------------
Tápanyagtartalom
Fehérje  : 12.0 g
Zsír     : 6.0 g
Rost     : 0.3 g
Só       : 0.3 g

Most ugyanez absztrakt modellként:

In [79]:
# Létrehozok egy absztrakt modellt
model = pyo.AbstractModel()

# Deklarálom a változók indexelésére használt halmazokat
model.I  = pyo.Set() # Ingredients
model.Nl = pyo.Set() # Nutrients with lower bound
model.Nu = pyo.Set() # Nutrients with upper bound

# Deklarálom a célfüggvényt mint paraméter
model.c = pyo.Param(model.I)

# Deklarálom az együtthatómátrixot és korlátozóvektort a tápanyagokhoz
model.A = pyo.Param(model.I, model.Nl | model.Nu)
model.b = pyo.Param(model.Nl | model.Nu)

# Deklarálom a döntési változókat, amelyeket az I halmaz elemei fognak indexelni
model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals)

# Definiálom, hogyan kell viselkednie a költségfüggvénynek
def cost(m):
    return pyo.summation(m.c, m.x)

# Definiálom, hogyan kell viselkednie egy felső tápanyagkorlátnak
def nutrient_upper_bounds(m, n):
    return sum(m.A[i,n] * m.x[i] for i in m.I) <= m.b[n]

# Definiálom, hogyan kell viselkednie egy alsó tápanyagkorlátnak
def nutrient_lower_bounds(m, n):
    return sum(m.A[i,n] * m.x[i] for i in m.I) >= m.b[n]
    
# Definiálom, hogyan kell viselkednie a nettó tömegre vonatkozó korlátnak
def net_weight_rule(m):
    return sum(m.x[i] for i in m.I) == 100

# Beállítom a költségfüggvényt, használva a korábban definiált szabályt
model.cost = pyo.Objective(rule=cost)

# Beállítom az alsó és felső tápanyagkorlátokat az Nl és Nu halmazok minden elemére
model.nutrientLBConstraint = pyo.Constraint(model.Nl, rule=nutrient_lower_bounds)
model.nutrientUBConstraint = pyo.Constraint(model.Nu, rule=nutrient_upper_bounds)

# Beállítom a nettó tömegre vonatkozó korlátot
model.netWeightConstraint  = pyo.Constraint(rule=net_weight_rule)

Az absztrakt modellt példányosítani kell konkrét adatokkal, én a nutrients.dat fájlban tároltam el a tápérték információt. A .dat fileok szintaxisáról: https://pyomo.readthedocs.io/en/stable/working_abstractmodels/data/datfiles.html#

In [67]:
!type nutrients.dat
set I := CSIRKE MARHA BIRKA RIZS KORPA ZSELATIN;
set Nl := feherje zsir;
set Nu := rost so;

param c :=
    CSIRKE   0.013
    MARHA    0.008
    BIRKA    0.010
    RIZS     0.002
    KORPA    0.005
    ZSELATIN 0.001;

param A : feherje zsir rost so :=
    CSIRKE   0.100 0.080 0.001 0.002
    MARHA    0.200 0.100 0.005 0.005
    BIRKA    0.150 0.110 0.003 0.007
    RIZS     0.000 0.010 0.100 0.002
    KORPA    0.040 0.010 0.150 0.008
    ZSELATIN 0.000 0.000 0.000 0.000;

param b :=
    feherje  8
    zsir     6
    rost     2
    so       0.4;
In [80]:
# Készítek egy példány a modellből a "nutrients.dat" file adataival
instance = model.create_instance(filename="nutrients.dat")

# Inicializálok egy solvert, ami az Xpresst fogja használni
solver = pyo.SolverFactory("xpress")

# Megoldom a modellt
solver.solve(instance)

# Kiírom a költséget, összetételt és tápanyagtartalmat
print("Költség  :", instance.cost(), "$")
print("--------------------")
print("Összetétel")
for i in instance.I:
    print(i, ":", instance.x[i](), "g")
print("--------------------")
print("Tápanyagtartalom")
for j in instance.Nl:
    print(j, ":", instance.nutrientLBConstraint[j](), "g")
for j in instance.Nu:
    print(j, ":", instance.nutrientUBConstraint[j](), "g")
WARNING: Extracting subsets for Set b_index, which is a SetOperator other than
    a SetProduct.  Returning this set and not descending into the set
    operands.  To descend into this operator, specify
    'subsets(expand_all_set_operators=True)' or to suppress this warning,
    specify 'subsets(expand_all_set_operators=False)'
WARNING: Extracting subsets for Set b_index, which is a SetOperator other than
    a SetProduct.  Returning this set and not descending into the set
    operands.  To descend into this operator, specify
    'subsets(expand_all_set_operators=True)' or to suppress this warning,
    specify 'subsets(expand_all_set_operators=False)'
Költség  : 0.52 $
--------------------
Összetétel
CSIRKE : 0.0 g
MARHA : 60.0 g
BIRKA : 0.0 g
RIZS : 0.0 g
KORPA : 0.0 g
ZSELATIN : 40.0 g
--------------------
Tápanyagtartalom
feherje : 12.0 g
zsir : 6.0 g
rost : 0.3 g
so : 0.3 g

PuLP

A PyOMO-hoz hasonlóan a PuLP is egy modellezésre fejlesztett könyvtár, amely külső solvereket használ az optimalizálási feladat megoldására. Az egyik előnye, hogy a CBC (COIN-OR Branch and Cut) solver, ami egy open source solver, alapértelmezetten települ vele együtt, így nem kell azzal bajlódnunk, hogy telepítsünk további solvereket (de megtehetjük természetesen). Ehhez útmutató: https://coin-or.github.io/pulp/

Fontos különbség a Pyomo és a PuLP között, hogy a PuLP-ban minden modell konkrét, illetve nem támogatja nem-lineáris modellek leírását, míg a Pyomo igen.

A modell alapegysége itt az LpProblem, illetve az LpVariable (amire előírható, hogy egészértékű legyen). A modellhez a változókból képzett lineáris kifejezéseket lehet hozzáadni korlátként és célfüggvényként.

Nézzük meg a keverési feladatot példaként!

In [55]:
!pip install pulp
Requirement already satisfied: pulp in c:\users\peter dobrovoczki\anaconda3\envs\notebook\lib\site-packages (2.5.1)
In [56]:
# Definiálom az összetevőket
ingredients = ["CSIRKE", "MARHA", "BIRKA", "RIZS", "KORPA", "ZSELATIN"]

# Definiálom a tápanyagokat
nutrientsLower = ["fehérje", "zsír"]
nutrientsUpper = ["rost", "só"]

# Definiálom a költségeket
costs = {
    "CSIRKE"  : 0.013,
    "MARHA"   : 0.008,
    "BIRKA"   : 0.010,
    "RIZS"    : 0.002,
    "KORPA"   : 0.005,
    "ZSELATIN": 0.001,
}

nutrientPercentage = {}

# Definiálom a tápanyagtartalmakat
nutrientPercentage["fehérje"] = {
    "CSIRKE"  : 0.100,
    "MARHA"   : 0.200,
    "BIRKA"   : 0.150,
    "RIZS"    : 0.000,
    "KORPA"   : 0.040,
    "ZSELATIN": 0.000,
}

nutrientPercentage["zsír"] = {
    "CSIRKE"  : 0.080,
    "MARHA"   : 0.100,
    "BIRKA"   : 0.110,
    "RIZS"    : 0.010,
    "KORPA"   : 0.010,
    "ZSELATIN": 0.000,
}

nutrientPercentage["rost"] = {
    "CSIRKE"  : 0.001,
    "MARHA"   : 0.005,
    "BIRKA"   : 0.003,
    "RIZS"    : 0.100,
    "KORPA"   : 0.150,
    "ZSELATIN": 0.000,
}

nutrientPercentage["só"] = {
    "CSIRKE"  : 0.002,
    "MARHA"   : 0.005,
    "BIRKA"   : 0.007,
    "RIZS"    : 0.002,
    "KORPA"   : 0.008,
    "ZSELATIN": 0.000
}

# Definiálom a tápanyagkorlátokat
nutrientBounds = {
    "fehérje" : 8,
    "zsír"    : 6,
    "rost"    : 2,
    "só"      : 0.4

}
In [81]:
import pulp as pl

# Inicializálom a modellt, "KeverésiProbléma"-nak nevezem el
# és előírom hogy minimalizálni fogom a célfüggvényt
prob = pl.LpProblem(name="KeverésiProbléma",
                    sense=pl.LpMinimize)

# Definiálom a változókat, az indexeiket az 'ingredients' listából kapják,
# 0 alsó korláttal és mindegyik folytonos
x = pl.LpVariable.dicts(name="Összetevők",
                        indexs=ingredients,
                        lowBound=0,
                        cat="Continuous") # egészértékű változó esetén "Integer" lenne

# Költségfüggvény
prob += (pl.lpSum(costs[i] * x[i] for i in ingredients), "költség")

# Alsó korlátok a tápanyagokra
for n in nutrientsLower:
    prob += (pl.lpSum(nutrientPercentage[n][i] * x[i] for i in ingredients) >= nutrientBounds[n], n)

# Felső korlátok a tápanyagokra
for n in nutrientsUpper:
    prob += (pl.lpSum(nutrientPercentage[n][i] * x[i] for i in ingredients) <= nutrientBounds[n], n)
    
# Korlát a nettó tömegre
prob += (pl.lpSum(x[i] for i in ingredients) == 100, "nettó tömeg")

# Kiírom a jelenlegi LP-t
print(prob)

# Megoldom
prob.solve()

# Kiírom az optimális költséget, összetételt és tápanyagokat
print("Költség:", pl.value(prob.objective), "$")
print("--------------------")
print("Összetevők")
for v in prob.variables():
    print(v.name.split("_")[1], ":", v.varValue, "g")
print("--------------------")
print("Tápanyagok")    
for n in nutrientsLower:
    print(n, ":", pl.value(prob.constraints[n]) + nutrientBounds[n], "g")
for n in nutrientsUpper:    
    print(n, ":", pl.value(prob.constraints[n]) + nutrientBounds[n], "g")
KeverésiProbléma:
MINIMIZE
0.01*Összetevők_BIRKA + 0.013*Összetevők_CSIRKE + 0.005*Összetevők_KORPA + 0.008*Összetevők_MARHA + 0.002*Összetevők_RIZS + 0.001*Összetevők_ZSELATIN + 0.0
SUBJECT TO
fehérje: 0.15 Összetevők_BIRKA + 0.1 Összetevők_CSIRKE + 0.04 Összetevők_KORPA
 + 0.2 Összetevők_MARHA >= 8

zsír: 0.11 Összetevők_BIRKA + 0.08 Összetevők_CSIRKE + 0.01 Összetevők_KORPA
 + 0.1 Összetevők_MARHA + 0.01 Összetevők_RIZS >= 6

rost: 0.003 Összetevők_BIRKA + 0.001 Összetevők_CSIRKE + 0.15 Összetevők_KORPA
 + 0.005 Összetevők_MARHA + 0.1 Összetevők_RIZS <= 2

só: 0.007 Összetevők_BIRKA + 0.002 Összetevők_CSIRKE + 0.008 Összetevők_KORPA
 + 0.005 Összetevők_MARHA + 0.002 Összetevők_RIZS <= 0.4

nettó_tömeg: Összetevők_BIRKA + Összetevők_CSIRKE + Összetevők_KORPA
 + Összetevők_MARHA + Összetevők_RIZS + Összetevők_ZSELATIN = 100

VARIABLES
Összetevők_BIRKA Continuous
Összetevők_CSIRKE Continuous
Összetevők_KORPA Continuous
Összetevők_MARHA Continuous
Összetevők_RIZS Continuous
Összetevők_ZSELATIN Continuous

Költség: 0.52 $
--------------------
Összetevők
BIRKA : 0.0 g
CSIRKE : 0.0 g
KORPA : 0.0 g
MARHA : 60.0 g
RIZS : 0.0 g
ZSELATIN : 40.0 g
--------------------
Tápanyagok
fehérje : 12.0 g
zsír : 6.0 g
rost : 0.30000000000000004 g
só : 0.3 g